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ABSTRACT 

A thermodynamically-based work potential theory for modeling progres- 
sive damage and failure in fiber-reinforced laminates is presented. The cur- 
rent, multiple-internal state variable (ISV) formulation, referred to as enhanced 
Schapery theory (EST), utilizes separate IS Vs for modeling the effects of damage 
and failure. Consistent characteristic lengths are introduced into the formulation 
to govern the evolution of the failure ISVs. Using the stationarity of the total 
work potential with respect to each ISV, a set of thermodynamically consistent 
evolution equations for the ISVs are derived. The theory is implemented into 
a commercial finite element code. The model is verified against experimental 
results from two laminated, T800/3900-2 panels containing a central notch and 
different fiber-orientation stacking sequences. Global load versus displacement, 
global load versus local strain gage data, and macroscopic failure paths obtained 
from the models are compared against the experimental results. 


INTRODUCTION 

Development of accurate, predictive progressive damage and failure analysis 
(PDFA) tools is imperative for further advancement of composite structures. 
As the predictive capability of PDFA tools improves the potential design space 
for a composite component is expanded. Physical testing of materials and com- 
ponents will always be the most crucial aspect of design validation; however 
virtual testing can help focus, and reduce, the necessary testing on only the 
most promising designs. To meet the predictive requirements needed for virtual 
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testing, a PDFA tool must be numerical objective and must incorporate the 
physics of the underlying damage and failure mechanisms in the composite. 

A thermodynamically-based, work potential theory, known as Schapery the- 
ory (ST), was developed for modeling matrix microdamage in fiber-reinforced 
laminates (FRLs)[l, 2, 3, 4], With ST, a single ISV is used to mark the degree of 
microdamage in the material. However, separate microdamage functions are em- 
ployed to degrade the transverse and shear moduli, yielding coupled anisotropic 
damage evolution. [5] and [6] extended the formulation to include the effects 
of transverse cracking by adding an additional ISVs and predicted the evolu- 
tion of microdamage and transverse cracking in coupon laminates analytically. 
[7] implemented this extended formulation in a numerical setting to simulate 
the failure of a buffer strip- reinforced, center- notched panel (CNP). However, 
due to the cumbersome nature of the evolution equations, the microdamage 
and transverse cracking evolution equations were decoupled to arrive at a more 
efficient implementation. Since no characteristic length is introduced into the 
formulation, the theory produces mesh-dependent results in a computational 
setting. 

The ST formulation is further extended here, resulting in the enhanced 
Schapery theory (EST), to include the effects of macroscopic transverse and 
shear matrix cracking, as well as fiber breakage, using an approach that differs 
from [5, 6, 7]. A deliberate distinction between damage and failure is made in the 
theory. The traditional ISV used in ST is maintained to model microdamage. 

Upon failure initiation, the element domain is no longer considered a con- 
tinuum, and a smeared crack approach is used to model the embedded discon- 
tinuities [8, 9]. Three new ISVs, which incorporate the characteristic length of 
the finite element and operate similar to traditional tensorial damage variables, 
dictate the directional evolution of the failure due to fiber breakage and trans- 
verse cracking mechanisms. Since EST is implemented within FEM, transverse 
cracking is considered a discrete phenomenon leading immediately to failure at 
the element scale (i.e., one macroscopic crack per element). This does not, how- 
ever, preclude the progressive evolution of cracks at the global scale that has 
been demonstrated to occur in multi-angle laminates [10]. The EST formula- 
tion presented herein offers non-linear progressive damage coupled with mesh 
objective, post-peak strain softening. EST incorporates physics-based methods 
for controlling the damage evolution within the material by using the measured 
relationships between transverse/shear degradation and the dissipated potential 
resulting from matrix microdamage. Moreover, failure evolution is governed by 
the material fracture toughness associated with fiber breakage and transverse 
cracking. 

Finally, EST is verified against experimental results for two center-notched 
panels (CNPs). Global load versus deflection data, local strain gage data, as 
well as observed failure mechanisms obtained from experiments performed at the 
NASA Langley Research Center (LaRC) and reported in [11, 12] are compared 
to numerical results. 
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Figure 1. Typical stress-strain curve, containing pre-peak nonlinearity and post-peak strain 
softening, showing the total elastic (We) and total dissipated ( Ws ) potentials. 


ENHANCED SCHAPERY THEORY 

The previously developed ST ([3, 6, 7, 13, 14, 15]) is extended to accom- 
modate mesh objective, post peak strain softening. Separate ISVs are used to 
govern the evolution of matrix microdamage, transverse (mode I) matrix fail- 
ure, shear (mode II) matrix failure, and fiber breakage (mode I). The first and 
second laws of thermodynamics are enforced, establishing thermodynamically 
consistent evolution laws for progressive matrix microdamage, as well as post- 
peak failure. The following sections give a brief overview of the formulation of 
this work potential theory. For complete details, please refer to [16, 17, 18] 


Thermodynamically-based Work Potential Framework 

As a material is loaded, a measure of the work potential facilitates modeling 
structural changes in the material, such as microcracking, which affect the elas- 
tic properties of the material. Energy that is not dissipated is recovered when 
the structure is unloaded, and the magnitude of energy recovered is contin- 
gent upon the degraded, elastic properties at the previously attained maximum 
strain state. It is assumed, upon subsequent reloading, that the material be- 
haves linearly, exhibiting the elastic properties observed during unloading, until 
the material reaches the preceding maximum strain state. After this state is 
achieved, structural changes resume, affecting degradation of the instantaneous 
elastic moduli of the material. This process is shown in the uniaxial stress-strain 
curve displayed in Figure 1. The shaded area above the unloading line represents 
total dissipated potential Ws, and the triangular area underneath is the total 
elastic strain energy density We- It is assumed that the material behaves as a 
secant material and there is no permanent deformation upon unloading. This 
a reasonable assumption for FRLs [5]; however, plastic deformation can also be 
incorporated, if necessary [3]. Extension to treat viscoelastic and viscoplastic 
response is outlined in [19]. 

Both We and Ws are functions of a set of ISVs, S m , (m = 1, 2, M). These 
ISVs account for any inelastic structural changes in the material. It is shown in 


[2, 3] that the total work potential is stationary with respect to each ISV. 


dW T 

~dS^ 


( 1 ) 


Additionally, [20] utilized the second law of thermodynamics to establish the 
inequality: 

Sm > 0 (2) 

where the overdot represents a time derivative and suggests that “healing” is 
not allowed for a material undergoing structural changes. Equations (1) and (2) 
form the foundation of a thermodynamically-based work potential theory for 
modeling non-linear structural changes in a material exhibiting limited path- 
dependence. 


Multiple ISV Formulation of ST to Account for Multiple Damage 
and Failure Mechanisms 

Due to the generality of the evolution equations, Equations (1) and (2), 
the work potential theory can account for any number and type of structural 
changes that may occur in a material. This is especially useful for modeling 
progressive damage in composites because the heterogeneity of the composite, 
and multiaxiality of the local fields, enables multiple damage mechanisms to 
arise during a typical loading history. 

The present EST formulation assumes that three major intralaminar mech- 
anisms are responsible for all observed non-linearities in the stress-strain curve 
of a composite lamina: matrix microdamage, matrix macroscopic cracking, and 
axial fiber failure. Each of these mechanisms can be accommodated by partition- 
ing the total dissipated energy density, Ws, into portions associated with each 
mechanism. Although delamination is a prominent failure mechanism in lami- 
nated composites, theories and methods to address this mechanism have been 
developed to a fairly mature state [21, 22]. The authors have left the coupling 
of in-plane and out-of-plane damage mechanisms as a future exercise involving 
the placement of cohesive zone elements between plies, yet recent preliminary 
results have been promising [23]. 

It has been reported that, under tension at room temperature, thermoset 
epoxy matrix composites exhibit minimal plastic deformation or viscous effects [1, 
5]. Thus, matrix microdamage is the primary cause of observed non-linearity in 
the stress versus strain response of many polymer matrix composites (PMCs) 
up to localization of microdamage into more severe failure mechanisms, such 
as transverse cracking, fiber breakage, kink band formation, or delamination. 
Microdamage can be considered as the combination of matrix microcracking, 
micro- void growth, shear banding, and fiber-matrix debonding. Figure 2 shows a 
typical uniaxial response of a material exhibiting microdamage evolution, where 
the recoverable energy potential is given by W and the potential dissipated into 
evolving structural changes associated with microdamage is given by S. 
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Figure 2. Typical stress-strain curve with a positive-definite tangent stiffness exhibiting 
microdamage, showing the elastic {W) and irrecoverable ( S ) portions. 

Typically, matrix microdamage continues to grow until the onset of more 
catastrophic failure mechanisms initiate. It should be noted that we have ex- 
plicitly distinguished between damage and failure in the following manner: 

Damage - Structural changes in a material that manifest as pre-peak non- 
linearity in the stress-strain response of the material through the degra- 
dation of the secant moduli. 

Failure - Structural changes that result from damage localization in a material 
and manifest as post-peak strain softening in the stress-strain response of 
the material. 

It is noted that this demarcation is not inclusive nor unique with respect 
to other types of phenomena in a broad range of materials that can also affect 
the secant moduli. The purpose here is to clearly point out demarcations with 
respect to the pre-peak (positive tangent stiffness) and post-peak (non-positive 
tangent stiffness) regime observed in matrix dominated damage and failure as 
it applies to fiber reinforced polymer matrix composites. 

These demarcations of damage and failure differ slightly from those typically 
used in previous works. Often, damage would be defined as any phenomena that 
affected the secant stiffness of the material (regardless of whether it yielded a 
positive or negative tangent stiffness), and failure was considered the final point 
on a stress-strain curve, after which, the stress would drop suddenly to zero [24, 
25, 26, 11, 27, 28, 29, 30]. By partitioning damage and failure, as specified in this 
work, a characteristic length can be introduced into the formulation of the failure 
evolution, as necessary, while the damage evolution can remain independent of 
the size. 

Here, three major failure mechanisms, which are distinct from the microdam- 
age modes, are considered: transverse (mode I) matrix cracking, shear (mode 
II) matrix cracking, and axial (mode I) fiber fracture. These failure modes are 
consistent with the in-plane failure typically observed in PMC laminates. It is 
assumed that the evolution of these mechanisms yields an immediate reduction 
in the load-carrying capability of a local subvolume where the mechanism is ac- 
tive. Three ISVs are used to account for mode I matrix cracking, mode II matrix 


cracking, and mode I fiber failure, respectively: S™, S J}, and S{. These ISVs are 
defined completely in a later section, and are taken to be the potentials required 
to advance structural changes associated with these failure mechanisms. 

At any given state, the total dissipated energy density, IT's, can be calculated 
as a sum of energy dissipated through the aforementioned damage and failure 
mechanisms, given by the four ISVs. 

W s = S + S? F + S? IF + S{ (3) 

According to the first law of thermodynamics, the total work potential (ignoring 
thermal dissipation) is given by the sum of the elastic strain energy density and 
the potentials associated with each of the damage or failure mechanisms. 

W T = W E + S + S? + + S{ (4) 


where W E is the elastic strain energy density. Invoking the stationarity principle, 
Equation (1), 

dW E „ dW E 


dS 


= -1 


dS{ 


= -1 


( 5 ) 


dW E = OWe 
dS™ dS ™ 

and the Second Law of Thermodynamics, Equation (2), gives: 


S > 0 S f F > 0 

sy F > o st if > o 


( 6 ) 


Equations (5) and (6) constitute the evolution equations for damage and failure 
in a material associated with matrix microdamage, matrix cracking, and fiber 
breakage in tension. 


Use of Traction-Separation Relationships to Define the Failure Po- 
tentials 

In the previous publications [5, 6, 7], it was assumed that the transverse 
cracking affected the relationships between stress and strain. However, the exis- 
tence of a macroscopic crack invalidates the assumption of a continuum. Here, 
it is presumed that failure arises from the evolution of cohesive cracks within the 
continuum, and the ISVs associated with failure (axial, transverse, and shear) 
influence the relationship between traction on the crack faces and the crack-tip 
opening displacement. The satisfaction of a Hashin-Rotem (H-R) failure crite- 
rion [31] indicates the material behavior transitions from that of a damaging 
continuum to that of a cohesive crack, and the essential fields become traction 
and separation, rather than stress and strain (see Figure 3). 

Once a cohesive crack initiates in the continuum, opening of the crack yields 
a reduction in traction on the crack faces at the crack tip, [32], If subsequently 
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Figure 3. Schematic showing the transition form a continuum to a cohesive zone due to the 
initiation of macroscopic cracks. The essential, constitutive variables switch from stress and 

strain to traction and separation. 



(a) Mode I fracture. 



(b) Mode II fracture. 


Figure 4. Triangular traction versus separation which dictates the behavior of cohesive 
cracks embedded in the continuum. The total area under the traction-separation law 
represents the material fracture toughness lC . 


the crack is closed, it is assumed that traction at the crack tip will unload 
linearly towards the origin of the traction versus separation law (see Figure 4). 
The strain energy release rate (SERR) G J M is calculated as: 

= / tU< ( 7 ) 

where j indicates the material (fiber / or matrix m), M represents the corre- 
sponding mode (mode I or mode II), S J M is the crack tip opening displacement 
in mode M and material j, and t J M is the corresponding traction at the crack 
tip. 

It should be noted, the traction-separation laws exhibited in Figure 4 do not 
require any initial, fictitious, pre-peak stiffness because the cracks are embedded 
within a continuum. This is an advantage over the use of DCZM elements which 
do require an initial stiffness because these interfacial elements do not actually 
represent physical material within the model and must attempt to simulate 
initially perfect bonding between adjacent material domains [33, 34], If set 
incorrectly, these fictitious stiffnesses can cause numerical problems [35]. 

Using the traction-separation laws in Figure 4 and assuming the cracks are 
not closing, the dissipation potentials in an element resulting from macroscopic 


cracking are related to the SERRs using suitable element dimensions and re- 
moving the energy potential that may be recovered upon closing of the crack. 

G J +f XI rim j-mxm rim j.m rm 

= I t I°I cm _ r J_ _ l I °I cm Lr II _ t II°II /ON 

°I .(0+90°) 9/ (0+9O°) ’ AO) 9/ (0) ’ AO) 9/ (0) ’ ^ 

where, li e+9 ° ^ is the length of a line running perpendicular to fiber direction in 
the element that intersects two edges of the element and the integration point, 
and le is the length of a line that is parallel to the fiber direction in the element 
that intersects two edges of the element and the integration point. 


EST Evolution Equations for a Fiber-Reinforced Lamina 

To arrive at the evolution equations for the four ISVs, the elastic strain en- 
ergy density must be defined for a material which may contain cohesive cracks. 
Therefore, the elastic strain energy We is comprised of a contribution from the 
continuum W and any possible cohesive cracks W 3 M . The plane stress, elastic 
strain energy density in the continuum is defined as 

R = -(-Ene^ + E 22 (S)e 2 22 + Gi 2 (S) 7 i 2 2 ) + < 512^11622 (9) 

where stress in the laminae are related to strain assuming plane stress conditions 

[36]. 

Note that only the transverse and shear moduli (E 22 and G\ 2 ) are functions 
of S since matrix microdamage only accrues in the matrix of the laminae. The 
Poisson’s ratio is assumed to evolve such that the quantity Q \ 2 = E 22 Vi 2 remains 
constant; however, this restriction can be relaxed if deemed necessary. The 
degraded moduli are related to the virgin moduli (E220 and G+20) and the ISV 
through a set of microdamage functions (e s (S) and g s (S )) that are obtained 
from three uniaxial coupon tests [2, 5, 6]. 


E 22 — E 220 e s (S ) 


( 10 ) 


Gi 2 — G 120 g s (S) (11) 

Degrading E 22 and G ± 2 exclusively is consistent with the intralaminar damage 
typically observed in PMC laminates. 

The elastic strain energy density of the cohesive cracks are defined as the 
recoverable energy per unit crack surface area smeared over the entire element. 
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(12) 
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(14) 


Substituting Equations (9) - (11) into Equation (5) gives the evolution equa- 
tion for the microdamage ISV. 

I + 7?2Gi2 °l|) = ~ 3S * (15) 

The use of a reduced ISV S r = S i has been employed in Equation (15). [5] has 
shown that the use of this reduced ISV yields polynomial forms of the micro- 
damage functions in Equations (10) and (11). Thus, Equation (15) becomes a 
polynomial equation that can be readily solved for S r for a given strain state 
(^22, 7l2)- 

Once failure initiates, the effects of failure supersede the effects of micro- 
damage and evolution of S ceases. The triangular traction-separation laws (see 
Figure 4) can be employed to evaluate Equations (12)-(14), then they can be 
substituted into Equations (5), which can be solved for the degraded moduli of 
the lamina. 
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where E* n and G\ 2 are the degraded transverse and shear moduli, due to micro- 
damage, when the H-R criterion is satisfied. 

To transition from crack tip opening displacement to integration point strain 
(which is needed to perfrom the calculations at an integration point within a 
FEM model), it is assumed that following failure initiation the strains are related 
to the crack tip opening displacements by 


i<»+*n fll = ((»+*>”> f c + S f 

(19) 

if e 22 = if 4 + (’T 

(20) 

if 712 = If ' lg + 2 ST, 

(21) 


where e ( f \ , e^ 2 , and 7^ are the strains when the H-R failure criterion is satisfied. 
Equations (19)-(21) imply that the strain in the continuum remains at the values 
obtained when failure initiates, and that any incremental change in the global 
strain after failure initiation is used wholly to advance the crack tip opening 
displacement. 

The negative tangent stiffness of the stress-strain response necessary for post- 
peak strain softening to occur imposes a restriction on the maximum allowable 
element size, as shown by [8]. 
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The analyst must be careful to ensure the dimensions of any failing elements are 
smaller than the conditions given in Equations (22) and (23). 

In summary, the satisfaction of the H-R failure criterion marks the transition 
from evolving microdamage to failure to macroscopic cracking. Prior to failure 
initiation, Equation (15) is used to calculate the microdamage reduced ISV S r , 
and the failure ISVs Sf, S™, and Sf} remain zero. Equations (10) and (11) 
are used to calculate the degraded transverse and shear moduli. Subsequent to 
matrix failure initiation, microdamage growth is precluded, and S r remains at 
S*, the value of S r when the H-R criterion for the matrix was satisfied. The 
degeneration of the transverse and shear moduli, resulting from matrix trans- 
verse and shear cracking, is calculated using Equations (17) and (18). Finally if 
the H-R criterion for the fiber is satisfied, the axial modulus is calculated using 
Equations (16) as fiber breakage evolves in the element. Once the material mod- 
uli have been calculated using the appropriate evolution equations, the stresses 
can be updated accordingly. 


EXAMPLE - CENTER NOTCHED PANELS SUBJECTED TO 
UNIAXIAL TENSION 

Experimental Details 

Two center-notched panel (CNP) configurations were tested at the NASA 
Langley Research Center (LaRC) [11, 12]. The geometrical details of the panel 
and testing boundary conditions are presented in Figure 5. End tabs affixed to 
the panel were clamped and a vertical, tensile displacement (in the ^/-direction) 
was applied to the top tab using a servo- hydraulic testing machine. The bottom 
tab was fixed preventing any ^-displacement of the bottom boundary of the gage 
section. The gripped tabs also prevented any displacement in the rc-direction at 
the top and bottom boundaries of the gage section. 

The panels were comprised of laminated T800/3900-2 carbon fiber/toughened 
epoxy composites. Three different lay-up configurations were tested; however 
one of the configurations exhibited significant delamination. Since the focus of 


Figure 5. 





Geometry, boundary conditions, and strain gage (Sg) locations of CNPs tested at 

NASA LaRC [11]. 


TABLE I. T800/3900-2 LAY-UP CONFIGURATIONS USED IN CNP TESTS AT NASA 

LaRC. 


ID 

Stacking Sequence 

Thickness (in.) 

Laminate- 1 

[0°] 12 

0.078 

Laminate-2 

[45°/0%45 o /0 o /90°],s 

0.065 


this work is modeling in-plane damage and failure mechanisms, this configura- 
tion is not considered here. The two remaining configurations are presented in 
Table 1. The first lay-up, Laminate-1, consists of twelve 0° plies, and the second, 
Laminate-2, is a symmetric, multi-angle lay-up with 40% of |45°|, 40% of 0°, and 
20% of 90° laminae. 

Several strain gages where affixed to the test panel, labeled Sg-1 through 
Sg-4 in Figure 5. Sg-1 was placed in the center of the panel, 1.5” above the 
notch. Sg-2 was placed 1.5” above the notch tip. Sg-3 was attached in front of 
the notch, 0.5” from the free edge, and Sg-4 placed at the notch tip. Global load 
versus displacement data, and local strain gage data were reported in [11, 12], 
along with a post-test C-Scan of Laminate- 1 and photograph of Laminate- 1 and 
photograph of Laminate-2. 

Finite Element Model Details 

EST was integrated into the Abaqus FEM software using the user material 
Fortran subroutine UMAT[37]. The linear elastic properties of T800/3900-2 used 
in the FEM models were taken from [11], The shear microdamage function g s 
utilized in Equation (11) was obtained from the shear stress-shear strain (T 12 - 712 ) 
response of [45°/-45 0 ] 3 s angle-ply T800/3900-2 coupon tests as recommended 
by [5]. The transverse, tensile and compressive microdamage functions were 
inferred by scaling the coefficients of the microdamage curves presented by [5] for 
AS4/3502 by the ratio of the virgin transverse modulus of T800/3900-2 to that 
of AS4/3502, as the stress-strain curves of the T800/3900-2 coupon laminates 
necessary to characterize e s were not available. The microdamage curves are 
presented in Figure 6. 



(a) Shear microdamage function 
obtained from [±45°]s angle-ply 
laminate. 



(b) Transverse tension and compression 
microdamage functions obtained by 
scaling data for AS4/3502 in Ref. [5]. 


Figure 6. Microdamage functions for T800/3900-2 used in FEM models. 
TABLE II. FAILURE PARAMETERS FOR T800/3900-2. 


Property 

Value 

Property 

Value 

X T 

0.021 

nf 

Wc 

1026 

in. 

Y t 

0.0092 

U IC 

2.39 

in. 

Z 

0.0075 

^ IIC 

6.78 

in. 


The axial mode I, transverse mode I, and shear mode II critical cohesive 
strains, and fracture toughness’ are given in Table II. The matrix mode II co- 
hesive critical strains ( Z ) and the fracture toughnesses (Gt} c ) were calibrated 
using data from Laminate-1. The critical strains were estimated using the ratio 
of the matrix strengths of IM7/8552, [38] and the calibrated critical shear strain 
in Table II. It was assumed that the characteristic material length of the mode 
I transverse crack band mechanism of T800/3900-2 and IM7/8552 were compa- 
rable. Thus, the ratio of the mode I matrix fracture toughness to the strength 
squared was multiplied by the assumed critical mode I matrix strain and the 
transverse modulus to arrive at the value of G ™ c used in the simulations and 
tabulated in Table II. Laminate- 1 did not exhibit any axial failure; so, the fiber 
mode I parameters {Xt and G IC ) were calibrated such that the ultimate load 
from the simulation of Laminate- 2 corresponded with the ultimate load reported 
by [11] for Laminate-2. 


Results - Laminate-1 

Global load P versus displacement A of a 4” section of Laminate- 1 is com- 
pared to results from the EST simulation in Figure 7. Very good agreement 
between the model and the experimental results is achieved. The response of 
the specimen appears to be linear until near 8,000 lbf., where the specimen be- 
gins deforming non-linearly. The EST simulation captures the initiation and 
progression of the global nonlinearity accurately. This panel was not loaded 


until catastrophic failure; hence, the data presented in Figure 7 represents load 
versus displacement data prior to the ultimate load of the specimen. 



Figure 7. Applied load versus displacement of a 4” section for Laminate-1. 

Local strain gage data (global load P versus local ^/-direction strain e yy ) from 
Laminate- 1 is plotted with the results from the EST FEM model in Figure 8; 
refer to Figure 5 for locations of strain gages. Strain relaxation is observed in 
the gage farthest away from the notch: Sg-1, shown in Figure 8a. The mode 
I and mode II matrix failure parameters in EST were calibrated such that the 
model demonstrates the same transition into strain relaxation at this location 
and at a similar global applied load. This load, taken as the splitting load, is 
8,250 lbf. in the experiment and 8,210 lbf. in the model. The transition to strain 
relaxation is more abrupt in the experiment as evidenced by the sharp knee in 
the load-strain curve, whereas, the transition in the model is more gradual. 

The data from the experiment and simulation for Sg-2, which is located 1.5” 
directly above the notch, are presented in Figure 8b. The model predicts less 
strain at Sg-2, for a given load, than the experiment, but the non-linear trends 
are very similar. This gage lies directly in front of the splitting crack path, 
shown in Figure 5, and it is not realistic to expect perfect agreement in areas 
experiencing high levels of damage and failure because of idealizations used to 
model the evolution of cracks in the simulations. 

Figure 8c displays data for Sg-3, located in front of the notch near the free 
edge. Very good agreement between the experimental and simulation results 
are exhibited. The model accurately captures the non-linear evolution of strain, 
away from the highly damaged regions, as a function of applied load. 

Finally, results for Sg-4 (located directly at the notch tip) are given in Figure 
8d and includes the experiment and simulation display of axial strain relaxation. 
As with Sg-2, Sg-4 shows less strain for a given applied load. However, the load 
at which the strain at Sg-4 relaxes in both the experiment and model correlate 
well, in accordance with the splitting load. Again, the relaxation response of 
the experiment is discontinuous, but the model exhibits continuous behavior. 

A C-Scan of the failed Laminate- 1 specimen is displayed in Figure 9. Four 
splitting cracks can be observed, shown in green, propagating outward from the 
notch tip, parallel to the loading direction, towards the gripped edges. Contour 
plots of the normalized microdamage obtained from the simulation are presented 






(c) Sg-3. (d) Sg-4. 


Figure 8. Applied load versus local strain for Laminate-1. 


in Figure 10a at 16,400 lbf. In these plots, S r is normalized by the maximum 
achievable value. The regions of maximum damage are localized to small regions, 
along the same crack path observed in Figure 9, embedded in a more widespread 
domain and exhibiting less severe microdamage. Figure 10b shows the shear 
failure degradation factor D at 16,400 lbf. The shear failure localizes into 
crack bands that are a single element wide and progress equivalently to the 
splitting cracks observed in the experiment. 


Results - Laminate-2 

Numerical results for applied load versus displacement of a 4” section of 
Laminate-2 are presented in Figure 11. The experimental ultimate load 15,300 
lbf. correlates well (axial failure parameters were calibrated to obtain an ulti- 
mate load that most closely matched the experimental data) with the ultimate 
load obtained from the model, also 15,300 lbf.. The global response up to failure 
is nearly linear and failure occurs suddenly and catastrophically. 

Figure 12 compares the applied load versus strain gage results from the model 
to the data from the experiment. Sg-1 and Sg-2 exhibited similar behavior; 
the strain increases until the ultimate load is obtained, after which the strain 
relaxes abruptly. The experimental data and numerical results both display 



Damage Path (Splitting) 


Figure 9. C-Scan of failed Laminate-1 specimen [11]. 
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Figure 10. Normalized matrix microdamage contour s ^ r ax and Matrix shear failure D™j in 

Laminate- 1 at P = 16,400 lbf.. 


this behavior. The model exhibits slightly more strain, for a given load, prior 
to ultimate failure. At Sg-3, the model predicts strain localization after the 
ultimate load is achieved. The gage data shows a slight reduction in strain 
as the load drops; however, the gage was placed directly in the crack path 
and may have been damaged when the panel failed. The model results and 
experimental data for Sg-4 exhibit similar trends, but the strain gage shows a 
large degree of nonlinearity at the notch tip. This observed nonlinearity was 
attributed to local interlaminar stresses near the notch free edge which caused 
some local delaminations [11]. Since the focus of this work was modeling in- 
plane damage mechanisms, these effects are not captured; however, the model 
could be easily extended to incorporate delamination by placing DCZM elements 
between continuum shell layers [12], 

A photograph taken of the failed, Laminate- 2 specimen is presented in Figure 
13. The photograph shows that two macroscopic cracks initially propagate from 
the notch tip towards the free edges, perpendicular to the applied load, in a 
self-similar fashion. Eventually, the cracks turn and proceed towards the free 
edge at an angle. [11] claim, supported by visual image correlation displacement 
data, that there was some eccentricity in the specimen alignment, which resulted 
in deviation from self-similar crack growth. 

Contours representing the normalized microdamage in the 45° layer is pre- 
sented in Figure 14 after the panel has completely failed and lost all of its load 
carrying capability. Although further matrix microdamage evolution is pro- 
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Figure 11. Applied load versus displacement of a 4” section for Laminate-2. 

hibited in elements that have failed (transverse/shear or axial), in the other 
elements that have not failed, matrix microdamage evolution continues. Nearly 
the entire 45° and -45° layers reach a normalized microdamage level of 0.18. The 
0° and 90° (not shown) plies exhibitted similar microdamage patterns; however, 
low levels of microdamage were more widespread in the 90° ply. 

Figure 15a shows the fiber failure path once the specimen has completely 
failed. All of the layers, except the 90° layer (not shown), show self similar 
cracks propagating from the notch tips towards the free edges of the panel. The 
angled crack path shown in Figure 13 was not reproduced because the eccentric 
loading (suspected in the test) was not introduced into the simulation; therefore, 
the crack growth remained self-similar. 

A high degree of transverse matrix failure can be seen in the axial crack path 
in the 45° ply in Figure 15b. This is also observed in the other plies, except the 
0° layer. Finally, Dy} is presented after the specimen has failed in Figure 15d. 
A highly degraded region has localized in the axial crack path. 

CONCLUSION 

A thermodynamically-based, work potential theory for damage and failure 
in composite materials, referred to as enhanced Schapery theory (EST), was 
developed. A marked distinction between damage and failure was introduced. 
Damage was considered to be the evolution of mechanisms that cause structural 
changes in the material such that the non-linear tangent stiffness tensor remains 
positive definite. Failure was taken to be the effect of structural changes in the 
material that result in loss of positive definiteness of the tangent stiffness matrix 
and post-peak strain softening. 

Separate internal state variables (ISVs) were used to account for damage 
and three in-plane failure mechanisms. The relationship between the transverse 
and shear moduli of the lamina were related to the microdamage ISV through 
a pair of experimentally-obtainable microdamage functions. Evolution of the 
failure ISVs is based upon traction-separation laws (which are a functions of 
the appropriate fracture toughnesses) and a characteristic element length. No 
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Figure 12. Applied load versus local strain for Laminate-2. 


mixed-mode law is used to couple mode I and mode II failure, but one can be 
readily implemented. Typically, the existence of a non-positive definite stiffness 
tensor would result in pathologically mesh dependent solutions; however, in 
EST, mesh objectivity is ensured by incorporating a characteristic length scale 
into the failure evolution. 

Two center- notched panels composed of T800/3900-2 were tested under ten- 
sile loading at NASA LaRC. Global load versus displacement and global load 
versus local strain gage strain data were compared to results obtained from 
FEM models utilizing EST. A single data point from each laminate was used to 
calibrate the EST failure parameters, and very good correlation was achieved 
for both laminates. The remaining strain gage and/or load-deflection results 
obtained from the EST simulations represent predictions and matched well with 
experimental data. Furthermore, damage and failure paths predicted by the 
models corresponded with the experimental results for both laminates exam- 
ined. This lends confidence to the approach presented in this paper for predict- 
ing acreage damage and failure of laminated composite structures. 



Figure 13. Photograph of failed Laminate-2 specimen [11]. 
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Figure 14. Normalized matrix microdamage contour in 45° layer of Laminate-2 after 
specimen has lost load carrying capability. 
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Figure 15. Failure in selected plies Laminate-2 after specimen has lost load carrying 

capability. 
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